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ABSTRACT 


V 

The  search  for  improved  performance  of  U.  S.  Navy  ships 
has  led  to  more  complex  propulsion  systems  consisting  of 
multiple,  interacting  inputs.  Classical  control  theory  does 
not  effectively  exploit  these  interactions.  Modern  control 
theory  provides  a  systematic  method  of  dealing  with  multiple 
interacting  inputs  to.  achieve  improved  system  performance. 
One  of  the  the  most  highly  developed  modern  control  tech¬ 
niques  is  the  linear  quadratic  regulator  ( LQR)  method. 
Essential  to  the  application  of  this  method  is  the  formula¬ 
tion  of  a  state  space  description  of  the  plant.  In  this 
paper  a  nonlinear  dynamic  propulsion  system  model  is  devel¬ 
oped  from  experimental  data  and  used  to  formulate  a  state 
space  model.  (  7/trsr-s.J 
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The  search  for  improved  performance  of  U.  S.  Navy  ships 
has  led  to  increasingly  complex  marine  propulsion  systems. 
Controllable  inputs  to  this  system  now  include  fuel  flow 
rate,  engine  "inlet  guide  vane  and  stator  vane  position, 
bleed  air  selection,  and  propeller  pitch  angle.  Current 
control  strategies  applied  to  these  propulsion  plants,  and 
classical  control  techniques  in  general,  do  not  take  into 
account  the  interaction  between  these  inputs.  In  contrast, 
modern  control  techniques  (MCT)  provide  a  systematic  method 
to  achieve  improved  system  performance  when  dealing  with 
multiple,  interacting  inputs.  Specifically,  modern  control 
theory  methods  provide  the  following  benefits  not  found  in 
classical  control  methods  applied  to  multiple  input, 
multiple  output  (MIMO)  systems: 

(1)  Effective  treatment  of  coupled  input  interactions  to 
improve  performance, 

(2)  Rigorous  treatment  of  stability  questions, 

(3)  Systematic  control  design  which  reduce  iteration  and 
the  need  for  extensive  intuition  and  experience  in 
the  control  design  process. 

The  most  extensive  application  of  modern  control  theory 
to  date  is  the  F100  Turbofan  Multivariable  Control  Synthesis 
Program,  sponsored  by  the  Air  Force  Aero  Propulsion  Lab  and 
Nasa  Lewis  Research  Center.  [Ref.  1:  p.  43]  .  The  results  of 
this  program  demonstrated  that  modern  control  theory  tech¬ 
niques  provide  an  orderly,  effective,  systematic  approach  to 
controller  design  for  multiple,  interacting  input  systems. 

Current  work  at  the  Naval  Postgraduate  School,  Monterey, 
is  aimed  at  investigating  the  application  of  modern  control 
techniques  to  U.  S.  Navy  ship  propulsion  plants.  The 
initial  phase  of  this  effort  involves  the  application  of 
modern  control  techniques  to  a  low  power  gas  turbine  test 


facility  located  at  the  school.  In  the  context  of  this 
effort  the  goals  of  this  thesis  were: 

(1)  Develooment  of  an  accurate  digital  computer  model  of 
the  propulsion  test  facility  which  includes  all 
significant  plant  nonlinearities  and  dynamic  effects. 

(2)  Develop  a  linear  (state-space)  model  of  the  propul¬ 
sion  plant. 

The  nonlinear  model  is  a  valuable  tool  in  controller 
design.  First,  it  provides  a  means  to  test  control  strat¬ 
egies  without  risking  damage  to  the  actual  plant.  Second, 
it  provides  a  cost  effective  alternative  to  extensive 
controller  tests.  Finally,  in  this  work  the  nonlinear  model 
provided  the  basis  from  which  a  linear  (state  space)  model 
was  derived. 

The  state  space  model  of  the  propulsion  plant  is  essen¬ 
tial  for  future  controller  design  work  at  NPS  using  the 
linear  quadratic  regulator  (LQR)  technique,  which  is  the 
most  highly  developed  modern  control  method.  [Ref.  2: 
p. 653] . 
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This  thesis  is  organized  into  eight  chapters.  In  the 
following  chapter  a  description  of  the  test  facility  at  the 
Naval  Postgraduate  School  (NPS)  is  given.  Also  in  that 
chapter,  a  conceptual  model  of  the  plant  is  developed.  The 
plant  is  first  divided  into  functional  components.  The 
significant  plant  dynamics  are  identified,  and  the  number  of 
plant  components  is  ’reduced  so  that  only  the  degree  of 
complexity  necessary-  to  represent  the  significant  plant 
processes  is  retained.  Finally,  the  component  interactions 
are  defined  by  identifying  component  inputs  and  outputs. 
The  coiflponent  interactions  link  the  components  together  and 
form  a  basis  for  quantitative  modeling  of  the  plant. 

In  Chapter  4  quantitative  (  versus  conceptual  )  compo¬ 
nent  modeling  is  described.  Using  experimental  data,  the 
input/output  relations  for  each  component  ( as  defined  in  the 
previous  chapter  )  are  constructed  in  equation  form. 

In  Chapter  5  the  individual  component  equations  are 
joined  together  to  form  a  steady  state  plant  model. 

In  Chapter  6  the  differential  equations  governing  the 
plant  dynamics  (  identified  in  Chapter  3  )are  introduced 
into  the  steady  state  model.  In  this  way  a  nonlinear 
dynamic  plant  model  is  developed. 

In  Chapter  7  the  state  space  model  is  derived  from  the 
nonlinear  dynamic  model. 

Chapter  8  contains  conclusions  and  recommendations  for 
further  work. 


A.  PLANT  DESCRIPTION 

The  test  facility  (  hereafter  referred  to  as  'the 
plant' )  consists  of  a  Boeing  model  502-6A  175  horsepower  gas 
turbine  engine  and  a  Clayton  17-300  water  dynamometer. 
Figure  3.1  is  a  schematic  of  the  test  facility. 

The  gas  turbine  can  be  subdivided  into  the  gas  generator 
section  and  the  power  'take-off  section.  The  gas  generator 
section  consists  of  a  single  stage  centrifugal  compressor,  a 
dual  can  combustor,  an  accessory  drive  section,  and  a 
single-stage  axial  flow  high  pressure  turbine  (HPT).  The 
power  take-off  section  consists  of  a  single-stage  axial  flow 
free  power  turbine  (FPT).  Fuel  enters  the  combustor  via  the 
fuel  control,  which  consists  of  a  flyball  governor  and  an 
acceleration  limiter.  The  fuel  control  setting  is  adjusted 
by  an  electric-motor  driven  control  lever.  The  hot  gases 
drive  the  high  pressure  turbine.  The  high  pressure  turbine, 
in  turn,  drives  the  compressor  and  accessory  drive  section. 
The  accessory  drive  section  extracts  power  from  the  high 
pressure  turbine  to  drive  the  fuel  pump  and  governor,  the 
oil  pump,  and  a  tachometer  generator.  The  free  power 
turbine  extracts  energy  from  the  hot  gas  stream  and  drives 
the  dynamometer.  There  is  no  mechanical  connection  between 
the  high  pressure  turbine  and  the  free  power  turbine. 

The  water  dynamometer  acts  as  a  power  absorption  unit. 
The  power  turbine  torque  and  speed  are  adjusted  by  varying 
the  amount  of  water  in  the  dynamometer.  This  is  analogous 
to  changing  the-  pitch  on  a  controllable  pitch  propeller. 
Water  enters  and  leaves  the  dynamometer  via  electric  motor 
driven  load  and  unload  valves. 


B.  THE  CONCEPTUAL  PLANT  MODEL 

In  formulating  a  conceptual  plant  model  the  following 
issues  must  be  addressed: 

•  definition  of  plant  boundaries, 

•  identification  of  significant  plant  dynamics, 

•  extent  to  which  the  plant  must  be  divided  into  compo¬ 
nents  and  how  the  components  selected,  and 

•  definition  of  component  inputs  and  outputs. 


These  issues  are  discussed  below. 

1.  Plant  Boundaries 

Selection  of  plant  boundaries  is  important  since 
this  selection  determines  the  plant  inputs.  Because  the 
long  range  goal  of  this  project  is  the  implementation  of  an 
improved  control  system,  the  existing  fuel  control  was 
excluded  from  the  plant  model.  The  plant  boundary  was 
located  downstream  of  the  existing  fuel  control,  and  actual 
fuel  flow  to  the  combustor  (versus  fuel  command  voltage)  was 
established'  as  one  plant  input. 

A  model  of  the  dynamometer  was  developed  in  earlier 
work  by  Johnson  (Ref.  3].  This  model  accurately  describes 
the  behavior  of  the  dynamometer  during  loading  conditions 
(water  addition  to  the  dynamometer),  but  is  less  accurate 
during  unloading  conditions.  The  source  of  the  dynamometer 
model  inaccuracies  is  thought  to  involve  the  unload  valve 
behavior.  In  order  to  avoid  introduction  of  the  unload 
valve  inaccuracies  into  the  propulsion  model  developed  in 
this  study,  the  load  and  unload  valves  were  placed  outside 
the  plant  boundaries.  Thus,  actual  water  to  and  from  the 
dynamometer  (versus  water  command  voltage)  was  established 
as  the  second  plant  input. 

2.  Plant  Components 

Breaking  the  plant  into  components  facilitates  the 
identification  of  causual  relationships  and  plant  dynamics. 
Further,  the  plant  components  provide  the  foundation  of 
plant  model  development. 


The  selection  of  components  is  a  cut  and  try 
process.  A  first  cut  at  component  identification  is  shown 
in  Figure  3.  2. 

This  selection,  was  based  on  a  functional,  basis  and  assump¬ 
tions  about  the  significant  processes  occuring  within  the 
plant.  The  minimum  number  of  components  necessary  to 
account  for  these  processes  is  sought.  If  the  resulting 
model  is  insufficiently  accurate,  then  some  significant 
process  has  either  been  overlooked  or  improperly  described. 
In  either  case,  the  selection  of  components  must  be 
reevaluated. 


Paramount  to  model  accuracy  is  the  identification  of 
significant  plant  dynamics.  Fluid  momentum  and  compressi¬ 
bility,  heat  transfer,  energy  storage,  rotor  inertia,  and 
combustion  effects  are  all  possibilities.  However,  dynamic 
effects  can  only  be  considered  significant  in  the  practical 
sense  if  their  time  constants  are  neither  much  shorter  nor 
extremely  longer  than  those  for  the  controller  actuators  and 
sensors.  Previous  work  by  Szuch  (Ref.  4:  p.  243]  in  the  area 
of  aircraft  gas  turbine  controls  indicated  that  fluid 
momentum,  compressibility,  and  energy  storage  dynamics  occur 
too  rapidly  to  be  controlled,  while  heat  transfer  dynamics 
occur  too  slowly  to  be  important  in  the  control  problem. 

The  importance  of  combustion  dynamics  deserves 
special  discussion.  Szuch  [Ref.  4:  p.  243]  and  DeHoff 
[  Ref.  5:  p.  274]  concluded  that  combustion  dynamics  were  of 
too  high  frequency  to  be  important  in  controls  considera¬ 
tions.  In  contrast,  Rubis  [Ref.  6:  p.  56]  discusses  signifi¬ 
cant  transient  effects  associated  with  engine  torque 
development  in  response  to  fuel  flowrate  changes.  These 
effects  could  be  due  in  part  to  combustion  related  delays. 
In  the  present  work  combustion  effects  were  initially 
neglected.  The  resulting  model  produced  excessive 
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Figure  3.2  Propulsion  Plant  Components 


accelerations  when  compared  to  experimental  data.  When 
combustion  effects  were  modeled  as  a  fuel  energy  lag 
accurate  results  were  achieved.  This  fuel  energy  lag 
represents  the  delay  between  the  time  when  the  chemical 
energy  in  the  fuel  passes  through  the  fuel  nozzles  and  the 
time  when  the  mechanical  energy  is  realized  at  the  high 
pressure  turbine. 

In  addition  to  the  fuel  energy  lag, 
the  most  significant  plant  dynamics  are  the  rotor  inertia 
effects.  Thus,  the  equations  governing  the  plant  dynamic 
behavior  are: 

NG  =  ( QH-QC-QA-QFRG)/JG 
NS  =  ( QF-QD-QFRD)/JD 
E/Mf  =  l/(tS  +  1) 

Where  NG  =  gas  generator  acceleration, 

NS  =  free  power  turbine  acceleration, 

-JG  *  gas  generator  inertia, 

JD  *  combined  free  power  turbine  and 
dynamometer  inertia, 

QH  -  high  pressure  turbine  torque, 

QC  =  compressor  torque, 

QA  =  accessory  drive  torque, 

QFRG  =  gas  generator  frictional  torque, 
QF  =  free  power  turbine  torque, 

QD  —  dynamometer  torque, 

QFRD  =  combined  free  power  turbine  and 

dynamometer  frictional  torque, 

Mf  =  measured  fuel  flowrate  at  the 
fuel  nozzles, 

E  =  mechanical  energy  applied  at  the 
high  pressure  turbine, 
t  =  time  constant  associated  with  the 
combustion  process. 


In  this  study  the  auxiliary  torque  (QD)  and  gas  generator 
frictional  torque  ( QFRG )  were  lumped  into  the  compressor 
torque  (QD).  Also,  the  dynamometer/power  turbine  frictional 
torque  (QFRD)  was  lumped  into  the  dynamometer  torque  (QD). 
Having  made  these  simplifications  the  governing  dynamic 
equations  become: 


NG  =  ( QH-QC ) / JG 


NS  =  (QF-QD)/JD 


Mf/E  =  l/(tS  +  1) 


(3.  1) 


(3.2) 


(3.3) 


4.  Model  Simplification 

Once  the  significant  dynamic  effects  are  determined 
and  the  governing  equations  identified,  the  plant  model  may 
be  simplified  by  reducing  the  number  of  components.  As 
noted  in  the  overview,  it  will  be  necessary  during  the  model 
development  to  define  the  input/output  relationships  for 
each  component  in  equation  form.  If  the  dynamic  effects  can 
be  lumped  into  isolated  components,  then  the  input/output 
equations  for  these  components  are  already  known,-  they  are 
the  governing  differential  equations  for  the  plant  dynamics. 
Further,  if  the  dynamic  effects  are  lumped  into  isolated 
components,  then  the  input/output  equations  for  the 
remaining  components  can  be  obtained  from  steady  state  data 
since  these  components  are  assumed  to  contain  no  dynamic 
effects.  This  approach  was  applied  to  the  current  modeling 
problem  and  the  reduced  component  model  of  Figure  3.3  was 
produced. 
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In  this  work  two  dynamic  effects  are  the  accelera¬ 
tions  of  the  gas  generator  and  the  power  turbine/ 
dynamometer.  By  lumping  the  entire  gas  generator  inertia 
into  the  gas  generator  shaft,  the  gas  generator  acceleration 
dynamics  were  isolated  to  that  component.  Similarly,  the 
power  turbine/dynamometer  dynamics  were  lumped  into  the 
power  turbine/dynamometer  shaft.  The  third  dynamic  effect, 
the  fuel  combustion  dynamics,  was  lumped  into  a  single 
component,  a  first  order  lag  between  the  fuel  flowrate  input 
and  the  high  pressure  .turbine. 

Since  no  dynamic  effects  were  considered  to  occur  in 
the  inlet  bell  or  exhaust  duct,  these  components  were 
combined  with  the  compressor  and  free  power  turbine, 
respectively. 

5.  Component  Input3 /Outputs 

The  next  step  in  formulating  the  conceptual  model  is 
determining  the  inputs  and  outputs  of  each  component. 
Multiport  analysis  using  signal  pairs  at  "ports"  of  power 
transfer  is  one  useful  method  for  studying  component 
interaction. 

At  the  mechanical -rotational  port  (gas  generator  and 
power  turbine/dynamometer  shafts  )  the  signal  pair  is 
torque,  Q,  and  rotational  speed,  N,  as  shown  in  Figure  3. 4. 
More  difficult  to  represent  is  the  thermofluid  power 
transfer  between  the  compressor,  high  pressure  turbine,  and 
free  power  turbine.  If  the  fluid  flow  were  incompressible 
the  port  would  be  represented  by  a  pres sure- volume  flowrate 
signal  pair  as  shown  in  Figure  3.  5. 

However,  this  simple  representation  is  not  adequate 
for  the  case  of  compressible  flow  since  it  does  not  account 
for  thermal  energy  transfer  via  the  fluid  internal  energy. 
One  could  represent  this  port  as  shown  in  Figure  3.6 
where  U  =  M  *  u,  with  M  the  mass  flowrate  and  u  the  fluid 
specific  internal  energy. 
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Figure  3.  4  Mechanical-Rotational  Port 


Figure  3. 5  Incompressible  Fluid  Port 

Representing  internal  energy  as  the  product  of 
specific  heat  ( Cv)  and  temperature,  the  port  could  be 
rewritten  as  shown  in  Figure  3.7. 

One  might  now  attempt  to  combine  effects  and  describe  the 
complete  thermofluid  power  transfer  as  shown  in  Figure  3. 8. 

However,  since  the  variables  P,V,M,  and  T  are 
related  through  the  equations  of  state,  it  is  redundant  to 
measure  all  four  in  the  modeling  process. 


Figure  3. 6  Thermal  Port 


Figure  3.  7  Alternate  Thermal  Port 

In  considering  which  variable  to  eliminate  it  should 
be  noted  that  since  the  direction  of  mass  flow  and  volume 
flow  must  be  the  same,  and  since  temperature  is  "carried 
along"  with  the  mass  flow,  these  three  variables  have  the 
same  signal  flow  direction.  Thus,  elimination  of  the  pres¬ 
sure  variable  was  considered  unwise  since  this  would  essen¬ 
tially  eliminate  the  two-way  component  interaction. 
Elimination  of  temperature  from  the  set  of  independent 


M 


P 


Figure  3.8  Combined  Port 

variables  was  considered  a  poor  choice  due  to  its  ease  of 
measurement  relative  to  either  mass  or  volume  flowrate. 
Finally,  since  mass  flowrate  is  conserved  while  volume  flow- 
rate  is  not,  and  since  this  conservation  might  lead  to 
simplification  in  future  measurements  and  calculations,  it 
was  decided  that  volume  flowrate,  V,  would  be  eliminated. 
The  resulting  power  port  is  shown  in  Figure  3.9. 
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Using  this  concept  of  the  thermodynamic  power  port 
the  complete  multiport  diagram  was  constructed  as  shown  in 
Figure  3.  10.  Note  that  variable  ambient  conditions  were 
eliminated  as  system  inputs  by  using  corrected  variables,  a 
commonly  employed  technique  in  gas  turbine  analysis. 
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IV. 


From  the  multiport  diagram  of  the  previous  chapter  the 
component  input/output  relationships  are  identified.  The 
next  step  in  the  modeling  process  is  to  obtain  quantitative 
expressions  for  these  relationships.  In  this  model  we  have 
assumed  that  the  gas  generator  and  power  turbine/dynamometer 
shafts,  and  the  fuel  energy  lag  contain  the  plant  dynamic 
effects.  The  quantitative  relation  for  these  components  is 
the  governing  differential  equations  3.  1  ,  3.2  and  3.3. 
Further,  since  it  has  been  assumed  that  the  remaining  compo¬ 
nents  contain  no  significant  dynamic  effect,  the  input/ 
output  equations  for  these  components  can  be  obtained  from 
steady  state  data.  The  following  describes  how  these 
equations  were  obtained. 


A.  DATA  ACQUISITION 

The  gas  turbine  and  dynamometer  is  instrumented  as  indi¬ 
cated  in  Table  1.  The  data  acquisition  system  is  controlled 
by  an  HP-85  personal  computer.  Temperature,  speeds,  and 
torques  are  taken  using  an  HP-3497A  Data  Acquisition/Control 
Unit.  Pressure  readings  are  taken  using  a  Pressure  Systems 
DPT-6400.  Fuel  flowrate  is  obtained  by  the  operator  from 
two  rotometers  and  entered  interactively  or  from  a  turbine 
flowmeter.  An  HP-82902M  flexible  disk  drive  provides 
program/data  storage  capability.  A  Digital  DECWRITER  IV 
printer  provides  hard  copy  output. 

Data  was  taken  at  93  operating  points  as  indicated  in 
Table  2.  These  points  were  selected  to  provide  full 
coverage  of  the  plant  operating  envelope.  At  each  operating 
point  temperature,  speed,  and  torque  values  were  sampled  30 
times  and  averaged.  Pressure  values  were  sampled  8  times 
and  averaged.  A  copy  of  the  acquisition’program  is  included 
as  Appendix  A. 


TABLE  1 

GAS  TURBINE/DYNAMOMETER  INSTRUMENTATION 
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B.  DATA  REDUCTION 


Data  reduction  took  place  in  two  phases.  In  the  first 
phase  average  values,  generator  torque,  air  mass  flowrate, 
and  corrected  values  were  computed.  These  calculations  are 
part  of  the  data  acquisition  program  (Appendix  A). 

In  the  second  phase,  curve  fits  were  obtained  for  the 
input/output  relations  of  each  turbine  component  using  the 
method  of  least  squares  [Ref.  7:p.  153]  The  program  which 

performed  the  least  squares  fit  is  included  as  Appendix  B. 

Using  this  program,  three  types  of  curve  fit  were  obtained 
for  each  input/output  relation.  The  first  is  a  "complete 
quadratic"  curve  fit  shown  in  equation  4.1, 

Y  =  C1*X12  +  C2*X1*X2  +  C3*X22  +  C4*X1  +  C5*X2  +  C6  (4.1) 

where  Y  =  output  (dependent)  variable, 

XI, X2  =  input  (independent)  variables, 

C1,C2,C3,C4,C5,C6,  =  constant  coefficients. 

The  second  type  of  curve  fit  obtained  was  a  "reduced 
quadratic",  so  called  because  the  cross  product  terms  (  ie. , 

X1*X2)  found  in  the  "complete  quadratic"  curve  fit  were 
excluded.  Equation  4.  2  is  an  example  of  this  format. 

Y  =  C1*X12  +  C2  *X22  +  C3*X1  +  C4*X2  +  C5  (4.2) 

The  third  curve  fit  type  was  a  linear  curve  fit. 

Equation  4.3  is  an  example  of  this  format. 

Y  =  C1*X1  +  C2*X2  +  C3  (4.3) 

In  each  case  the  result  of  the  curve  fit  program  is  the  jj 

coefficients  of  the  curve  fit  equation.  Because  the  magni¬ 
tude  of  some  variables  is  much  larger  than  others  ( ie. ,  NG  = 

30,000  rpm,  P4  =  17  psia  )  it  was  necessary  to  scale  the  S 
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TABLE  2 

GAS  TURBINE  DATA  ACQUISITION  SCHEDULE 


GAS  GENERATOR 
SPEED,  NG 
(RPM) 

19078. 73 

21128. 36 
.  22547.57 

23891.  30 
25280. 03 
26665.  44 
28726.  66 

30067. 85 
31525.  20 
33488.  27 

19254. 72 
21854.  58 
23270. 02 

24516. 74 

25983. 85 

27296. 90 

28796. 48 
30080. 00 

31543.90 
33490. 05 
35372. 62 

-  21112.06 

22577. 27 
23896. 71 

25259. 17 

26508. 37 

27917. 37 
29370. 06 
30694.  21 
32464.  30 

23925. 74 
25269. 44 

26637. 27 
27807. 57 
29131.  64 

30443. 17 
32379. 32 
34523.  04 
27344. 96 
28696. 43 

30057. 51 
31489.  43 
33358. 05 

30060. 52 
31528.  78 
33291. 47 


DYNAMOMETER 
SPEED.  NS 
(  RPM) 

549. 12 

557. 66 
562. 21 

566.  05 
565. 37 

567.  74 

576. 29 
580. 45 

583. 29 
593. 19 
946. 33 
954.  31 
957.  89 
954.  13 
968. 10 
964.  39 
971. 97 
989. 08 
977. 28 
981.  S6 
996.  06 

1441.88  . 
14-44.77 
1445. 22 

1442. 62 

1455. 67 
1452. 72 
1468. 75 
1460. 35 

1459. 62 
1934.  28 

1943. 66 
1944. 69 
1944.  34 
1933- 95 
1939. 51 
1974.  51 
1951. 02 

2432. 68 
2441.  53 
2436.  86 
2431.  06 
2440. 31 
2918. 94 
2927. 09 
2930. 93 


GAS  GENERATOR 
SPEED,  NG 
(  RPM) 

19400. 72 
21829. 68 
23208. 31 
24584.  86 
25922. 19 
27331.  84 

29382. 51 

30825. 73 
32533. 96 
19080. 00 
21160. 08 

22570. 12 
23968. 66 

25285. 25 
26635. 89 

28080. 53 
29448.  34 
30824.  83 
32476.  51 
34504.  89 
19088. 39 

21852. 49 
23202.79 
24476.83 

25919. 51 
27302. 78 
28612. 59 

30036. 55 
31568. 65 

33358. 12 
24574.  83 
25845.  15 

27287. 50 
28665. 43 

29985. 26 

31485. 56 
33144.  56 
35366. 35 
28000. 63 
29346. 11 
30736.  87 
32441.  26 
29370. 30  s 

30653. 77 
■  32413.38 

34366. 77 


DYNAMOMETER 
SPEED,  NS 
(RPM) 

563. 30 

557. 98 
563. 58 
569. 89 

573. 93 

569. 30 

580. 98 
587. 61 
584.  86 

935. 15 
958. 80 
952. 79 
953. 97 
959.  41 

965. 15 
969.07 
976.  55 
979. 37 

978. 31 

983. 94 
1428. 13 
1431.  24 
1436.  46 
1438. 68 
1443. 97 
1453. 57 
1446. 09 
1458. 07 
1456. 94 
1461. 41 
1953. 73 
1943. 32 
1946. 16 
1947. 56 
1954.  42 
1961.  59 
1938. 01 
1963. 19 
2421.  92 
2433. 45 
2429.  97 
2441.  16 
2922. 37 
2906.  41 
2917. 06 
2917.  24 


variables  to  prevent  algorithmic  singularities  during  the 
least  squares  solution.  This  was  accomplished  by  dividing 
each  variable  by  a  scaling  factor.  The  scaling  factors, 
shown  in  table  3,  were  selected  so  that  each  variable  had  a 
range  of  0.  0  to  1.0  when  scaled. 


TABLE  3 

SCALING  FACTORS 

Variable 

Scaling  Factor 

Gas  Generator  Speed,  NG 

36,000  rpm. 

Compressor  Torque,  Qc 

130  ft.  lb. 

Air  Mass  Flowrate,  Ma 

13,000  lb/hr. 

Compressor  Discharge 

Temperature,  T2 

800  deg.  R. 

Compressor  Discharge 

Pressure,  P2 

43.  0  psia. 

Fuel  Mass  Flowrate,  Mf 

240  lb/hr. 

• 

High  Pressure  Turbine 

Torque,  QH 

130  ft.  lb. 

Combined  Air/Fuel  Mass 
■Flowrate,  Maf 

13,000  lb/hr. 

High  Pressure  Turbine 

Discharge  Temp. ,  T4 

1800  deg.  R. 

High  Pressure  Turbine 

Discharge  Press. ,  P4 

20.  0  psia. 

Dynamometer  Speed,  NS 

3 , 000  rpm. 

Free  Power  Turbine 

Torque,  QF 

480  ft.  lb. 

A.  STEADY  STATE  MODEL  ALGORITHM 

With  the  input/output  component  relations  defined  in 

equation  form  the  next  step  was  to  link  these  relations 
together  to  form  a  steady  state  plant  model.  Figure  5. 1 
presents  the  flowchart  describing  the  steady  plant  state 
model  algorithm.  In  this  algorithm  the  program  user  inputs 
the  gas  generator  and /dynamometer  speeds  at  which  the  plant 
parameters  are  to  be  evaluated.  The  program  makes  an 
initial  guess  at  the  steady  state  fuel  flowrate,  Mf.  A 
guess  is  also  made  for  the  compressor  and  high  pressure 
turbine  discharge  pressures,  P2  and  P4.  The  program  uses 
these  assumed  values  to  calculate  the  compressor  and  high 
pressure  turbine  outputs.  The  computed  compressor  discharge 
pressure  is  compared  with  the  assumed  value.  If  the  differ¬ 
ence  between  assumed  and  computed  value  ceeds  the  specified 
tolerance,  the  value  of  P2  is  updated.  Otherwise,  the  power 
turbine  outputs  are  calculated  and  the  computed  high  pres¬ 
sure  turbine  discharge  pressure  is  compared  with  the  assumed 
value.  Convergence  within  the  specified  tolerance  is  again 
required.  If  this  check  is 

met  the  compressor  and  high  pressure  turbine  torques  are 
compared. 

These  torques  should  be  equal  in  steady  state.  If  they  are 
not  the  assumed  value  of  fuel  flowrate  is  updated  and  the 
entire  process  is  repeated.  The  steady  state  computer  model 
is  included  as  Appendix  C. 

B.  STEADY  STATE  MODEL  RESULTS 

The  steady  state  computer  program  was  tested  at  various 
points  in  the  plant  operating  envelope.  The  output  of  the 
computer  model  was  compared  with  the  raw  data  for  the  same 


operating  points.  The  results  showed  excellent  agreement 
between  the  steady  state  model  and  the  raw  data  at  most 
operating  points.  The  only  exception  to  this  was  at 
extremely  high  gas  generator  speeds.  A  typical  comparison 
is  shown  in  Table  4. 

As  indicated  by  the  flowchart  of  Figure  5.  1,  at  any 
given  gas  generator  speed  (NG)  and  dynamometer  speed  (NS) 
combination,  the  program  must  hunt  for  the  fuel  flow  rate 
that  will  produce  zero  torque  differential  between  the 
compressor  and  high  pressure  turbine.  An  investigation  into 
the  nature  of  the  torque  differential  as  a  function  of  fuel 
flowrate  (at  fixed  NG  and  NS)  provides  some  interesting 
insights  into  the  performance  characteristics  of  the  engine. 
As  an  example,  in  Figure  5.2  the  torque  differential  is 
plotted  versus  fuel  flow  for  NG  =  30,000  rpm  and  NS  =  1,100 
rpm. 

This  plot  has  several  interesting  features.  First,  note 
that  there  are  two  values  of  fuel  flowrate  that  lead  to  zero 
torque  differential.  This  *  suggests  that  at  a  given  NG/NS 
combination  there  are  actually  two  equilibrium  conditions. 
It  is  thought  that  the  upper  fuel  flowrate  represents  a 
state  of  inefficient  operation.  Whatever  the  source  of  the 
higher  fuel  flowrate  condition,  comparison  with  the  raw 
data  indicates  that  the  normal  operating  mode  of  the  turbine 
is  at  the  lower  fuel  flowrate. 

A  second  feature  of  this  curve  is  the  peak  torque 
differential.  This  peak  is  considered  to  be  a  significant 
characteristic  of  the  gas  generator.  It  is  thought  to  be 
indicative  of  the  maximum  driving  torque  difference  attai¬ 
nable  at  the  specified  NG/NS  combination.  Thus,  if  one 
wishes  to  accelerate  the  gas  generator  quickly,  this  plot 
indicates  the  limits  of  achieveable  acceleration  as  well  as 
the  optimal  fuel  input  to  achieve  the  greatest  acceleration. 
Clearly,  more  fuel  does  not  necessarily  lead  to  greater 


acceleration.  In  fact,  this  plot  suggests  that  too  large  a 
fuel  increase  can  lead  to  deceleration.  One  should  note, 
however,  that  this  condition  may  be  impossible  to  achieve 
practically,  since  it  presumes  that  this  large  fuel  change 
can  be  made  instantly,  without  changing  the  NG/NS  combina¬ 
tion  ( ie. ,  a  perfect  step).  The  fuel  energy  lag  dynamics 
seem  to  exclude  this  possibility  in  the  current  application. 
Further,  the  existing  fuel  control  devices  on  most  gas 
turbine  facilities  would  likely  prevent  this  condition  from 
being  observed. 
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The  next  step  was  to  develop  a  nonlinear  dynamic  model 
by  introducing  the  plant  dynamic  equations  into  the  steady 
state  model. 

In  order  to  implement  the  governing  dynamic  equations, 
the  gas  generator  and  power  turbine/dynamometer  inertias 
must  be  known.  Johnson  [Ref.  3:p.  56]  concluded  that  the 
combined  power  turbine'/dynamometer  inertia  (JD)  is  insensi¬ 
tive  to  dyno  water  weight  and  that  a  value  of  JD  =  0. 6738 
lb. ft. s2  was  valid  throughout  the  operating  range. 

The  following  section  describes  the  determination  of  the 
gas  generator  inertia.  Subsequently,  the  development  and 
results  of  the  nonlinear  dynamic  model  is  described. 

A.  GAS  GENERATOR  INERTIA 

The  technical  manual  for  the  Boeing  502-6A  gas  turbine 
engine  [Ref.  8: p.  6]  lists  the  gas  generator  inertia  (JG)  as 
0.11  in.  lb.  s2.  However,  it  was  unclear  whether  this  inertia 
value  included  the  accessory  gearbox  inertia.  Further,  the 
equipment  configuration  at  the  NPS  test  facility  is  somewhat 
different  than  the  standard  configuration  described  in  the 
technical  manual.  Because  the  gas  generator  inertia  is 
crucial  to  accurate  dynamic  performance  prediction,  it  was 
desirable  to  verify  the  technical  manual  value 
experimentally. 

In  order  to  experimentally  determine  the  gas  generator 
inertia,  the  gas  generator  inlet  bell  and  nose  cone  was 
removed.  A  lever  arm  was  attached  to  the  compressor 
impeller  using  existing  bolt  holes  intended  for  impeller 
removal.  [Ref.  8: p.  VII-10].  A  spring  was  attached  at  each 
end  of  the  lever  arm  and  secured  to  the  base  of  the  turbine. 
The  resulting  experimental  set-up  is  shown  in  Figure  6. 1. 


Finally,  a  potentiometer  was  attached  to  the  lever  arm  and 
aligned  with  the  impeller  shaft  centerline  to  permit 
measurement  of  impeller  angular  position. 

The  lever  was  deflected  and  released  from  rest.  Using  a 
strip  chart  recorder  the  oscillatory  motion  of  the  gas 
generator  was  recorded.  This  procedure  was  repeated  ten 
times  so  that  good  average  values  could  be  obtained. 

A  simplified  diagram  of  the  experimental  apparatus  is 
given  in  Figure  6.2,  where  K  is  the  effective  spring 
constant,  J  is  the  tptal  polar  mass  moment  of  inertia  about 
the  gas  generator  axis,  9  is  the  angle,  of  rotation,  and  c  is 
the  system  damping  due  to  friction.  From  this  simplified 
diagram  the  differential  equation  for  viscously  damped  free 
vibration  is  found  to  be: 


J9+  C  0  +  K  0  =  0 


(6.  1) 


The  solution  to  this  equation  [Ref.  9:p.  25-32]  for  the 
underdamped  case  reveals  the  frequency  of  damped  oscillation 
to  be: 


Wd  =  Wn  \jl-  c2 


(6.2) 


where  Wd  =  frequency  of  damped  oscillation 
Wn  =  natural  frequency 
C  =  damping  ratio. 

Using  average  values,  the  frequency  of  damped  oscilla¬ 
tion  was  determined  from  the  strip  chart  readings  to  be: 

Wd  =  67.  544  rad/sec 


Figure  6. 2  Simplified  Diagram  of  Experimental  Apparatus 

Also  from  the  solution  to  the  free  vibration  equation, 
the  following  relation  involving  the  damping  ratio  is  found: 

2  tt  z/\J  1-  S2  =  ( l/n)*ln(Xo/Xn)  (6.3) 

where  n  »  number  of  elapse’d  oscillatiions, 

Xo  =  original  oscillation  amplitude, 

Xn  =  amplitude  after  n  cycles. 

Again  using  average  experimental  values,  the  damping 
ratio  was  determined  from  equation  6. 3  to  be: 

C  =  0. 0325 

Using  the  experimentally  obtained  damping  ratio  and 
damped  frequency,  the  natural  frequency,  Wn,  was  determined 
from  equation  6.  2  to  be: 

Wn  =  67. 5799  rad/sec 

The  natural  frequency  of  the  system  is  also  given  by. 


For  this  system  the  effective  spring  constant,  K,  is: 


K  =  2*ks*R2  ( 6. 5) 

where  ks  =  individual  spring  constant, 

R  =  distance  from  gas  generator  axis  to 
spring  attachment  point. 

The  individual  spring  constants  were  experimentally 
determined  to  be  ks  =  6.891  lb. /in.  ,  with  R  =7.0  inches. 
Solving  for  the  effective  spring  constant  ,K,  using  equation 
6. 5  the  total  system  inertia,  J,  was  calculated  from 
equation  6.  4  to  be: 

J  =0. 14786  in. lb. s2 

The  total  inertia  is  equal  to  the  sum  of  the  individual 
inertia  effects  of  the  gas  generator  rotor  (JG),  lever  arm 
(JL),  and  springs  (JS): 

J  =  JG  +  JL  +  JS.  (  6.  6) 

The  lever  arm  inertia  was  calculated  as: 

JL  =  ml*l2/12  =  0.014457  in.  lb.  s2  (6.7) 

where  ml  =  mass  of  the  lever  arm, 

1  =  length  of  the  lever  arm. 

The  combined  inertia  of  both  springs  was  calculated 
using  Rayleigh'  s  method  to  be: 

JS  =  2*ms*R2/3  =  0.01914  in.  lb.  s2  (6.8) 

where  ms  =  mass  of  each  spring. 
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With  the  total  (J),  spring  (JS),  and  lever  (JL)  inertias 
determined,  the  gas  generator  inertia  was  found  using 
equation  6. 6  to  be: 

JG  =  0. 1143  in.  lb. s2 


B.  NONLINEAR  DYNAMIC  PROGRAM 

The  nonlinear  dynamic  program  was  formulated  using 
Discrete  Simulation  Language  (DSL).  The  flowchart 
describing  the  dynamic  program  algorithm  is  given  in  Figure 
6.3.  The  user  must  enter  the  fuel  flowrate  and  dynamometer 
water  weight  as  a  function  of  time  by  editing  the  program 
prior  to  execution.  Upon  execution  of  the  program  the  user 
interactively  enters  the  initial  gas  generator  and  dynamom¬ 
eter  speeds.  The  steady  state  program  is  then  called  to 
determine  the  equilibrium  value  of  fuel  flowrate,  Mfo,  and 
dynamometer  water  weight,  Wwo.  A  time  step  is  then  taken 
and  the  new  fuel  flowrate  and  water  weight  is  determined. 
The  dynamic  effect  of  the  fuel  energy  lag  is  then  computed. 
The  steady  state  program  is  again  used  to  determine  the 
compressor,  high  pressure  turbine,  free  power  turbine,  and 
dynamometer  torques  (QC,  QH,  QF,  QD  respectively).  These 
torques  are  entered  into  the  dynamic  equations  describing 
the  gas  generator  and  dynamometer  accelerations.  accelera¬ 
tions  are  integrated  to  obtain  speeds.  A  check  is  made  to 
determine  if  the  run  time  is  exceeded.  If  not,  time  is 
again  incremented  and  the  loop  repeats.  A  copy  of  the 
nonlinear  dynamic  program  is  included  as  Appendix  D. 

In  order  to  validate  the  nonlinear  dynamic  program  the 
propulsion  plant  test  facility  was  subjected  to  step  changes 
in  commanded  fuel  flowrate  voltage.  The  resulting  fuel 
flowrate,  gas  generator  speed,  and  dynamometer  speed  was 
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Figure  6.  3  Flowchart  for  Nonlinear  Dynamic  Program 


recorded  using  a  multichannel  strip  chart.  The  recorded  fuel 
flowrate  versus  time  was  entered  into  the  dynamic  program  in 
tabular  form  and  was  used  to  exercise  the  program.  Various 
acceleration  and  deceleration  tests  were  conducted.  The 
experimental  data  were  compared  with  the  output  of  the 
nonlinear  dynamic  program.  Figures  6.4,  6.5,  and  6.6  illus¬ 
trate  the  results  obtained.  In  this  case  an  increase  in 
fuel  flowrate  was  applied  to  accelerate  the  gas  generator 
from  25,000  to  29,500  rpm  and  the  dynamometer  from  960  to 
1090  rpm.  This  represents  a  large  transient,  covering 
nearly  one  third  of  the  gas  generator  operating  envelope. 

Figure  6. 4  shows  the  fuel  flowrate  transient  in  response 
to  a  step  change  in  commanded  fuel  flowrate  voltage.  The 
discontinuities  shown  result  from  the  discrete  sampling 
effects  of  the  digital  flowmeter,  and  the  unsteady  nature  of 
the  fuel  flow  at  the  location  of  the  flowmeter,  just  down¬ 
stream  of  the  fuel  control  valve.  Figures  6.  5  and  6.  6  show 
the  resulting  transients  of  the  gas  generator  and  dynamom¬ 
eter.  Experimental  data  is  plotted  along  with  the  nonlinear 
dynamic  model  results.  The  results  show  excellent  agreement 
between  the  data  and  the  model. 


NONLINEAR  MODEL  RESULTS 


NONLINEAR  MODEL  RESULTS 


Dynamometer  Response 


Among  modern  control  theory  techniques,  the  linear  quad¬ 
ratic  regulator  (LQR)  method  is  the  most  highly  developed. 
This  method  coordinates  multiple  inputs  simultaneously  and 
provides  a  straightforward  manner  in  which  the  feedback  gain 
matrix  can  be  manipulated  to  achieve  the  desired  system 
performance.  The  LQR  method  and  calls  for  the  system  to  be 
represented  in  state  space  form  as  shown  below: 

X  =  A*X  +  B*U  ( 7. 1) 

where  X  =  state  vector, 

U  =  input  vector, 

A  =  state  coefficient  matrix, 

B  =  input  coefficient  matrix. 

In  order-  to  arrive  at  the  state  space  representation, 
one  must  resort  to  perturbational  variables.  The  excursion 
of  any  variable,  X,  away  from  its  initial  condition  can  be 
represented  by: 

X  =  Xo  +  x  {1.2) 

where  Xo  =  the  initial  value, 

x  =  dX  =  the  perturbation  from 
the  initial  value, 

X  =  the  current  value. 

Any  variable  can  be  represented  in  this  manner. 

In  this  study  the  states  are  the  gas  generator  speed 
(NG)  ,  power  turbine/dynamometer  speed  (NS),  and  mechanical 
energy  resulting  fuel  combustion  (E).  This  selection  is 
mandated  by  the  dynamic  equations  of  3.1,  3.2,  and  3.3,  one 
state  per  derivative  term  (Ref.  10: p.  6651 .  The  plant 
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KV 


Si 


inputs  are  the  fuel  flowrate  (Mf)  and  the  dynamometer  water 
weight  (Ww).  Using  perturbational  variables  the  state  space 
equation  becomes: 

■';•)=  a /;:!+ b{::i  . 

What  remains  to  be  done  is  to  determine  the  elements  of 
the  'A1  and  'B'  matrices,  which  contain  the  coefficients  of 
the  state  equation  set.  We  can  write  these  elements  symbol¬ 
ically  as: 

all  =  drig/dng  al2  =  brig/ans  al3  =  dng/ae 

a21  =  iris/ang  a22  =  dris/ans  a23  =  aris/^e 

a31  =  de/bng  a32  =  fte/bns  a33  =  be/be 

bll  =  ing/bmf  bl2  =  brig/aww 

b21  =  bns/bmf  b22  =  bris/bww 

b31  =  ^e/amf  b32  =  oeA>ww  . 

In  order  to  arrive  at  these  coefficients  a  Taylor  series 
expansion  was  carried  out  on  each  component  input/output 
equation  retaining  only  first  order  terms.  The  results  is  a 
set  of  linear  equations  which  can  reduced  to  the  state  space 
form  given  above.  A  detailed  solution  for  these  coeffi¬ 
cients  is  included  as  Appendix  E.  It  is  important  to  note 
that  these  coefficients  vary  with  operating  point.  A 
subroutine  was  added  to  the  steady  state  program  which  eval¬ 
uates  these  analytic  coefficient  expressions  at  user  speci¬ 
fied  operating  points  (SUBROUTINE  PART  in  Appendix  C). 
Table  5  shows  how  the  'A'  and  'B’  matrices  vary  with 
operating  point. 

Comparison  between  the  state  space  model  and  nonlinear 
dynamic  model  was  conducted  by  subjecting  both  to  the  fuel 
flowrate  inputs  used  in  the  nonlinear  dynamic  program  valir 
dation.  Shown  in  Figure  7.  1  and  7.  2  is  the  response  of  the 
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nonlinear  and  state  space  models.  In  this  example  the  'A' 
and  ' B '  matrices  were  evaluated  at  the  initial  condition. 
As  expected  for  large  perturbations  such  as  this,  the  state 
space  model,  with  its  linear  assumptions,  does  not  accu¬ 
rately  describe  the  behavior  of  this  highly  nonlinear  plant. 
The  limitations  of  the  state  space  model  are  important  in 
that  they  indicate  the  limits  of  the  LQR  controller  design 
and 

help  to  define  necessary  transitions  between  linear  approx¬ 
imations.  When  an  accurate  global  dynamic  model  of  the 
plant  is  needed  for  control  testing,  the  nonlinear  model 
will  be  used. 


ison  of  State  Space  vs.  Nonlinear  Model 
Dynamometer  Response 


An  accurate  nonlinear  dynamic  computer  model  of  the 
Naval  Postgraduate  School  marine  propulsion  test  facility 
has  been  developed.  A  state  space  model  has  been  derived 
from  the  nonlinear  model. 

Prior  to  controller  design  the  fuel  flowrate  and  dyna¬ 
mometer  water  weight  actuators  must  be  accurately  modeled 
and  included  in  both  ’the  nonlinear  dynamic  program  and  the 
state  space  model.  In  conjunction  with  this  effort,  a  new 
fuel  control  valve  should  be  installed  which  will  allow  more 
direct  control  of  fuel  flowrate.  When  this  is  accomplished 
a  new  controller  should  be  designed  using  modern  control 
techniques,  specifically  the  LQR  method.  Performance  tests 
should  then  be  conducted  on  the  turbine  using  both  the 
present  (classical)  controller  and  the  modern  controller. 
These  tests  should  include  power  level  transients  as  well  as 
simulated  sea  state  conditions. 


DATA  ACQUISITION  PROGRAM 
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FOLLOWING  IS  A  INDEX  OF  THE  PARAMETERS  USED  IN  THIS  PROGRAM 
A  STATE  SPACE  COEFFICIENT  MATRIX. 
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THE  FOLLOWING  ENSURES  THAT  THE  OUTPUT  STAYS  IN  WITHIN  LIMITS. 
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CONSTRUCT  THE  COMPLETE  QUADRATIC  EQUATION. 
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THIS  SUBROUTINE  PRODUCES  OUTPUT  'T4'  FOR  THE  GIVEN  INPUTS. 


DIMENSION  X(5),C(21),Z(6),XR(5) 
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COEFFICIENTS  OF  THE  QUADRATIC  CURVE  FIT. 
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NONLINEAR  DYNAMIC  PROGRAM 
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APPENDIX  £ 

STATE  EQUATION  FORMULATION 


1.  Note  that  for  any  variable  Y, 

y  =  Yo  +  Ay 

where  Y.  =  final  value 

Yo  =  initial  value 
Ay  s  change  in  Y 

Assume  the  initial  condition,  Yo,  is  known.  Then  the  final 
state,  Y,  can  be  found  if  AY  is  known. 

Applying  a  Taylor  series  approximation  to  Y,  where  Y  is  a 
function  of  X  ( i. e.  Y  =  Y( X) )  yields, 

Y-Yo  =AY  =  Yo+Yo'(AX)+Yo"(AX)V2!  ♦.  .  .  Yon/n!  . 

For  very  small  changes  in  X,  a  linear  approximation  can  be 
made: 


Y  =  dY  =  Yo  +  Yo*dX. 

As  a  convenience  of  notation  the  following  substitution  is 
made: 


dY  =  y 

2.  If  Y  is  a  function  of  X  (  i.  e.  Y  =  Y(  X) ) ,  then  (making 
the  linear  assumption) 

JY/bX  =  oYo/oX  +  fty/bX. 

But  dYo/iX  =  0,  since  Yo  is  a  constant. 

Then  in  general, 

4Y/oX  =  by/ax. 
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3.  In  what  follows  the  state  space  equation  set  for  the 
propulsion  test  facility  is  derived  by  using  a  linear 
Taylor  series  approximation  for  each  of  the  input/output 
equations  used  in  the  nonlinear  model. 

NG  =  (QH-QC)/JG 

QH  =  QH(  Ma, T2 , E, NG, P4) 

QC  =  QC( P2,NG) 

rig  =  (3QH/JMa)ma  +  (iQH/iT2)t2  +  ( JQH/3E)e  + 

.  .  .*(  4QH/<»NG)ng  +  OQH/‘iP4)p4  - 
.  . . ( bQC/bP2)p2  -  (&QC/J>NG)ng 
rig  a  al*ma  +  a2*t2  +  a3*e  +  a4*p4  +  a5*ng  - 

. . . a6*p2  -  a7*ng  ( 1) 

Ma  »  Ma( P2 , NG) 

ma  =  ( ^ Ma/d P2  ) p2  +  (5Ma/&NG)ng 

=  bl*p2  +  b2*ng  (2) 

P2  a  P2(Ma,T2,E,P4,NG) 

p2  =  ( ^P2/bMa)ma  +  (hP2/J>T2)t2  +  ('iP2/iE)e  + 

. . . (iP2/^NG)ng  +  (&P2/bP4)p4 
p2  *  el*ma  +  e2*t2  +  e3*e  +  e4*p4  +  e5*ng  (3) 
T2  =  T2(  P2,NG) 

t2  =  (aT2/i)P2)p2  +  (dT2/0NG)ng 

=  cl*p2  +  c2*ng  (4) 

substitute  (2),  (4)  into  (3), 

p2  a  el( bl*p2+b2*ng)+e2(  cl*p2+c2*ng)+.  .  .  (5) 

e3*e+e4*p4*e5*ng 
P4  =  P4( Maf ,T4,NS) 

p4  =  ( fcP4/aMaf )maf  +  ($P4/'4T4)t4  +  ... 

(  6P4/'4NS)ns 

p4  =  dl*maf  ♦  d2*t4  +  d3*ns 
maf  =  ma  +  e 

p4  a  dl*( ma+e)  +  d2*t4  +  d3*ns  (6) 

T4  =  T4( Ma, T2 , E, P4, NG) 

t4  =  ( dT4/'4Ma)ma  ♦  (aT4/3T2)t2  ♦  (5T4/3E)e  ♦ 

.  .  .  (^T4/ANG)ng  ♦  (  oT4/i>P4)p4 
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(7) 


t4  =  fl*ma  +  £2*t2  +  f3*e  +  f4*p4  +  f5*ng 
substitute  (6)  into  (7)  and  solving  for  t4, 

t4(l-f4*d2)  =  ma*(  fl  +  f4*dl)  +  t2*£2  +  (8) 

.  .  .  e*{  t3+f4*dl )  +  f4*d3*ns  +  f5*ng 
let  gl  =  (  l-f4*d2 ) 

g2  =  { fl+f4*dl)/gl 
g3  =  f2/gl 
g4  =  (  f3+f4*dl)/gl 
g5  =  f 4*d3/gl 
g6  =  f5/gl 

then  t4  =  g2*ma  +  g3*t2  +  g4*e  +  g5*ns  +  g6*ng 

grouping  terms  in  (  5-) , 

p2( l-el*bl-e2*cl)  =  ng(  el*b2+e2*c2+e5 )  +.  .  . 
e*e3  +  e4*p4 

let  g7  »  ( l-el*bl-e2*cl) 

g8  =  ( el*b2+e2*c2+e5 )/g7 
g9  =  e3/g7 
glO  =  e4/g7 

then  p2  a  g8*ng  +  g9*e  +  gl0*p4  ( 10) 

substitute  (2),  (4)  into  (1), 

ng  =  al( bl*p2+b2*ng)  +  a2( cl*p2+c2*ng)  +  ... 

a3*e  +  a4*p4  +  a5*ng  -  a6*p2  -  a7*ng 

collecting  terms, 

rig  =  p2( al*bl+a2*cl-a6)  +  ng( al*b2+a2*c2  +  a5-a7 ) .  . 
. . . a3*e  +  a4*p4 

let  gll  =  ( al*bl+a2*cl-a6 ) 

gl2  =  ( al*b2+a2*c2+a5-a7) 
gl3  =  a3 
gl4  =  a4 

then  rig  =  gll*p2  +  gl2*ng  +  gl3*e+  gl4*p4  (11) 

substitute  (10)  into  (11)  and  collecting  terms, 

rig  =  ng( gll*g8+gl2 )  +  e*( gll*g9+gl3 )  +  ... 

p4*(  gll*gl0+gl4) 
gl5  =  (  gll*g8+gl2 ) 
gl6  =  (gll*g9+gl3) 


let 


gl7  =  ( gll*gl0+gl4) 

then  rig  =  gl5*ng  +  gl6*e  +  gl7*p4  (12) 

substitute  (9)  into  (6)  and  collect  terms, 

p4  =  ma(dl+d2*g2)  +  e(dl+d2*g4)  +  t2*d2*g3  +.  .  . 
ns(d2*g5+d3)  +  ng*d2*g6 


gl8  = 

( dl+d2*g2 ) 

gl9  * 

( dl+d2*g4) 

g20  = 

d2*g3 

g21  = 

( d2*g5+d3 ) 

g  22  = 

d2*g6 

then  p4  =  gl8*ma  +  gl9*e  +  g20*t2  +  g21*ns.  .  .  (13) 

+  g22*ng 

substitute  (2)  and  (4)  into  (13)  and  collect  terms, 

p4  =  p2( gl8*bl+g20*cl )  +  ng( gl8*b2+g20*c2+g22 ) .  .  . 
. .  .  +  gl9*e  +  g21*ns 
let  g23  =  ( gl8*bl+g20*cl ) 

g24  =  ( gl8*b2+g20*c2+g22 ) 

p4  =  g23*p2  +  g24*ng  +  gl9*e  +  g21*ns  ( 14) 

substitute  ( 10)  into  ( 14)  and  collect  terms, 

p4( l-g23*gl0)  *  ng(  g23*g8+g24)  +  e(  g23*g9+gl9 ) .  .  . 

. . .  +  g21*ns 

let  g25  =  (l-g23*gl0) 

g26  =  (  g23*g8+g24)/g25 
g27  =  ( g23*g9+gl9)/g25 
g28  =  g21/g25 

p4  =  g26*ng  +  g27*e  +  g28*ns  (15) 

substitute  (15)  into  (12)  and  collect  terms, 

rig  *  ( gl5+gl7*g26)  +  e(  gl6+gl7*g27 )  +  ns*gl7*g28 
let  g29  =  ( gl5+gl7*g26 ) 

g30  =  ( gl6+gl7*g27) 
g31  =  gl7*g28 

rig  =  g29*ng  +  g30*e  +  g31*ns  (16) 


NS  =  (QF-QD)/JD 
QH  =  QH( NS,Maf ,T4) 


QD  =  QD(NS,WW) 

ns  =  ( dQF/jNS)ns  +  ( }QF/3Maf )ma£  +  (;QF/jT4)t4 
. . . ( 3QD/3NS)ns  +  ( 5QD/jWW)ww 
ns  =  zl*ns  +  z2*maf  +  z3*t4  -  z4*ns  -  z5*ww 
raaf  =  ma  +  e 

ris  =  zl*ns  +  z2*ma  +  z2*e  +  z3*t4  ...  ( 

-  z4*ns  -  z5*ww 

substitute  (9)  into  (19)  and  collect  terms, 

ris  =  ns( zl+z3*g5-z4)  +  ma(z2+z3*g2)  +  ... 

e(z2+z3*g4)  +  t2*z3*g3  +  ... 
z3*g6*ng  -  z5*ww 
let  g32  =  ( zl+z3*g5-z4) 

g33  =  ( z2+ z3*g2) 
g34  =  ( z2+z3*g4) 
g35  =  t2*z3*g3 
g36  =  z3*g6 

ris  =  g32*ns  +  g33*ma  +  g34*e  +  .  .  .  ( 

g35*t2  +  g36*ng  +  g36*ng  -  z5*ww 
substitute  (2)  and  (4)  into  (21)  and  collect  terms, 
ris  =  ns*g23  +  p2(  g33*bl+g35*cl )  +  ... 

ng( g33*b2+g35*c2+g36)  + 
e*g34  -  z5*ww 

let  g37  =  ( g33*bl+g35*cl ) 

g38  =  ( g33*b2+g35*c2+g36 ) 

ris  =  g32*ns  +  g37*p2  +  g38*ng  +  .  .  .  ( 

g34*e  -  z5*ww 

substitute  (10)  into  (22)  and  collect  terms, 

ris  =  ns*g32  +  ng( g37*g8+g38)  +  e( g37*g9+g34) . . 
+  p4*g37*gl0  -  z5*ww 
let  g39  =■  (g37*g8+g38) 

g40  =  ( g37*g9+g34 ) 
g41  =  p4*g37*gl0 

ris  =  g32*ns  +  g39*ng  +  g40*e  +  . . .  ( 

g41*p4  -  z5*ww  * 

substitute  (15)  into  (23)  and  collect  terms. 


ns  =  ns(  g32  +  g41*g28 )  +  ng( g39+g41*g26 )  ... 

+  e( g40+g41*g27 )  -  z5*ww 
let  g39  =  ( g32+g41*g28 ) 

g40  =  ( g39+g41*g26 ) 
g41  =  ( g40+g41*g27) 

nS  =  g42*ns  +  g43*ng  +  g44*e  -  z5*ww 

E  =  MF/( tS  +  1) 

e  =  -e/t  +  mf/t  (...t  =  time  constant) 

4.  Equations  (16),  (24),  and  (25)  comprise  the  plant  state 
equations. 
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